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Abstract 

This paper presents a modified numerical scheme for a class of Frac- 
tional Optimal Control Problems (FOCPs) formulated in Agrawal (2004) 
where a Fractional Derivative (FD) is defined in the Riemann-Liouville 
sense. In this scheme, the entire time domain is divided into several sub- 
domains, and a fractional derivative (FDs) at a time node point is approx- 
imated using a modified Griinwald-Letnikov approach. For the first order 
derivative, the proposed modified Griinwald-Letnikov definition leads to a 
central difference scheme. When the approximations are substituted into 
the Fractional Optimal Control (FCO) equations, it leads to a set of alge- 
braic equations which are solved using a direct numerical technique. Two 
examples, one time-invariant and the other time- variant, are considered to 
study the performance of the numerical scheme. Results show that 1) as 
the order of the derivative approaches an integer value, these formulations 
lead to solutions for integer order system, and 2) as the sizes of the sub- 
domains are reduced, the solutions converge. It is hoped that the present 
scheme would lead to stable numerical methods for fractional differential 
equations and optimal control problems. 

Keywords: Fractional calculus, Riemann-Liouville fractional derivatives, mod- 
ified Griinwald-Letnikov approximation, fractional optimal control 

1 Introduction 

Optimal Control Problems (OCPs) appear in engineering, science, economics, 
and many other fields. An extensive body of work exists in the area of optimal 
control of integer order dynamic systems ( Hestenes (1966), Bryson and Ho(1975), 
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Gregory and Lin(1992)). It was shown recently that fractional derivatives pro- 
vide more accurate behavior of a dynamic system (see Podlubny (1999) and the 
references there in). Therefore, formulations and numerical schemes for optimal 
control problems which account for fractional dynamics of these systems would 
be necessary. In this work, we develop a modified numerical scheme for a class of 
Fractional Optimal Control Problems whose dynamics is described by Fractional 
Differential Equations. 

Agrawal (Agrawal (2004)) defines a Fractional Dynamic System (FDS) as a 
system whose dynamics is described by Fractional Differential Equations (FDEs) , 
and a Fractional Optimal Control Problem (FOCP) as an optimal control problem 
for an FDS. A general formulation for FOCPs was proposed in Agrawal (2004). 
As it can be seen from literature, there is no much work in the field of optimal 
control of FDSs. The formulations of FOCPs comes from Fractional Variational 
Calculus (FVC) which is an emerging branch of fractional calculus. 

Riewe (Riewe (1996), Riewe (1997)) was first to formulate a fractional vari- 
ational mechanics problem. Riewe's major focus was to develop Lagrangian 
and Hamiltonian mechanics for dissipative systems. Agrawal (Agrawal (2001)) 
presented an ad hoc approach to obtain the differential equations of fraction- 
ally damped systems. In Agrawal (2002), Agrawal presented fractional Euler- 
Lagrange equations for Fractional Variational Problems (EVPs). Klimek (Klimek 
(2001)) presented a fractional sequential mechanics model with symmetric frac- 
tional derivatives. In Klimek (2002), Klimek presented stationary conservation 
laws for fractional differential equations with variable coefficients. Dreisigmeyer 
and Young (2003) presented nonconservative Lagrangian mechanics using a gen- 
eralized function approach. In Dreisigmeyer and Young(2004) the authors show 
that obtaining differential equations for a nonconservative system using FDs may 
not be possible. 

The fractional Euler-Lagrange equation has recently been used by Baleanu 
and coworker to model fractional Lagrangian with linear velocities (Baleanu 
and Avkar(2004)), fractional metafluid dynamics (Baleanu (2004)), fractional 
Lagrangian and Hamiltonian formulations of discrete and continuous systems 
(Muslih and Baleanu (2005a), Muslih and Baleanu (2005b), Muslih et al.(2006), 
Baleanu and Muslih(2005)) and Hamiltonian analysis of irregular systems (Baleanu 
(2006)). Tarasov and Zaslavsky have used variational Euler-Lagrange equation to 
derive fractional generalization of the Ginzburg-Landau equation for fractal me- 
dia (Tarasov and Zaslavsky (2005)) and dynamic systems subjected to nonholo- 
nomic constraints (Tarasov and Zaslavsky(2006)). In (Agrawal (2004), Agrawal 

(2005) ), the fractional variational calculus is applied to deterministic and stochas- 
tic analysis of fractional optimal control problems. Rabei, Ajlouni and Ghassib 

(2006) develop suitable Lagrangian and Hamiltonian for a fractional dynamic 
system, which they transform to fractional Schrodinger's equation and solve it. 
Stanislavsky (2006) presents a Hamiltonian formulation of a dynamic system. 
Atanackovic and Stankovic (2007) present existence and uniqueness criteria for 
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problems resulting from fractional variational calculus. 

In this paper, we present a direct numerical scheme for a class of Fractional 
Optimal Control Problems (FOCPs) formulated in Agrawal (2004). The scheme 
uses a modified Griinwald-Letnikov definition to approximate a fractional deriva- 
tive. For a first order derivative, this approximation leads to a central difference 
formula. For simplicity in the discussion to follow, this formulation is briefiy 
presented here. Two examples are solved to demonstrate the performance of the 
algorithm. 

2 Fractional Optimal Control Formulation 

In this section, we briefly present a Hamiltonian formulation for an FOCP. Con- 
sider the following FOCP: Find the optimal control u{t) for a FDS that minimizes 
the performance index 

J(u)^ ['f{x,u,t)dt (1) 
Jo 

and satisfles the system dynamic constraints 

oD^x^g{x,u,t), (2) 

and the initial condition 

x{0) = xo, (3) 

where x{t) is the state variable, t represents the time, / and g are two arbitrary 
functions, and oDfx represents the left Riemann-Liouville derivative of order a 
of X with respect to t. For the definitions of fractional derivatives and some of 
their applications, see (Podlubny (1999), Magin (2006), and Kilbas, Srivastava 
and Trujillo (2006)). Note that the upper limit of the integration is taken as 
1. We consider < a < 1. Further, we consider that x{t), u{t), f{x,u,t) and 
g{x,u,t) are all scalar functions. These conditions arc made for simplicity. The 
same procedure could be followed if the upper limit of integration and a are 
greater than 1, and x{t), u{t), f{x,u,t) and g{x,u,t) are vector functions. 

It should be pointed out that in traditional integer-order optimal Control, 
Eq. (1) may also include terminal terms. Such terms lead to nonzero terminal 
condition at t = 1. For the FOCP considered here, our formulation would require 
fractional terminal terms, the meaning of which may not be clear. For this reason, 
the terminal terms are not included in Eq. (1). 

To find the optimal control we define a modified performance index as 

J(u) = C\H{x,u,t) -XoD^xUt, (4) 
Jo 

where H{x, u, A, t) is the Hamiltonian of the system defined as 

H{x, u, A, t) = f{x, u, t) + Xg{x, u, t), (5) 
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and A is the Lagrange multiplier. Taking variations of Eq. (jl]) and using ([5]), the 
necessary equations for the optimal control are given as 

.CfA = f. (6) 

and 

oDtx = ^ (8) 

Following the approach presented in Agrawal (2004), we also require that 

A(1)=0. (9) 

Equations ([6])- ([9]) represent the necessary conditions in terms of a Hamiltonian 
for the optimal control of the FOCP defined above. It could be verified that the 
total time derivative of the Hamiltonian as defined above is not zero along the 
optimum trajectory even when / and g do not explicitly depend on t. This is a 
departure from the integer order optimal control theory. 

In the discussion to follow, we shall strictly focus on the following quadratic 
performance index 

J(u) = - / [q{t)x'^{t) +r{t)u^]dt, (10) 
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where q{t) > and r{t) > 0, and the system whose dynamics is described by the 
following linear FDE, 

oDfx = a{t)x + h{t)u. (11) 

Using Eqs. ([6]) to ([H]), it can be demonstrated that the necessary Euler-Lagrange 
equations for this system are (see also Agrawal (2004)), 

oD'^x = a{t)x - r-\t)b'^{t)\, (12) 

tD'^\ = q{t)x + a{t)\, (13) 

and 

u = -r-\t)b{t)X. (14) 

Equations f|T2l) to f|T^ will be used to develop a direct numerical scheme for a 
FOCP. 

3 A Modified Numerical Scheme for FOCPs 

In this section, we define modified Griinwald-Letnikov approximations of frac- 
tional derivatives as 

oD-x(t,_i/2) = E z = 1, ■ ■ ■ , n. (15) 

^ 3=0 



n—i 

iDXt,+i/2) = E ^ = n - 1, n - 2, ■ ■ ■ , 0, (16) 

j=0 

where j = 0, 1, ■ ■ ■ , n are the coefficients. A recursive approach of computing 
ujj^^ is given as (Podlubny (1999)) 

/ - 1 , - (^ ^ + A , ,(°) , - 1 

'^0 — ~\ j I ^i-l' J — -'-5 

It can be shown that for a = 1, Eqs. f|T5l) and f|T6l) lead to 

(i3:(tj_i/2) _ Xj - Xj-i 
dt ~ h 

and 

dx{ti+l/2) _ Xj - Xj+i 

dt h 

which are essentially the central difference equations for the left and the right 
derivatives. 

To develop a numerical scheme, we divide the time domain [0, 1] into n equal 
parts, and approximate the fractional derivatives oD'^x and tD^X at the center 
of each part using Eqs. (ITSj) and (fT6l) . We further take x(ti_i/2) as an average 
of the two end values of the segment. Thus, a;(tj„i/2) = {xi-i + Xi)/2. We make 
similar approximations for A(ti_i/2), x{ti+i/2), and A(tj+i/2). Substituting these 
approximations into f|T2l) and f|T3l) . we obtain 

1 * 1 1 

I^E^j°^^»-i = -a{iih){xi^i +Xi) - -r'^{iih)b'^{iih){Xi^i + Xi) 
n ^.^Q Z Z 



^ = l,■■■,n (17) 

1 1 1 

T^^^j°'^-^i+j = :5'?(^2/i)(a;i+i + Xi) + -a{i2h){Xi-i + Xi) 
n z z 

i = n-l,---,0 (18) 

where ii = i — \ and 12 = i + \. Equations f|T7j) and f|T8l) represent a set of 2n 
linear equations in terms of 2n unknowns, which can be solved using a standard 
linear solver. One can also develop an iterative scheme in which one can march 
forward to compute Xi's and backward to compute Aj's to save storage space and 
perhaps computational time. 
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4 Numerical Examples 



To demonstrate the applicability of the formulation and to validate the numerical 
scheme, in this section we present numerical results for two problems, time invari- 
ant and time varying. For both problems, two types of studies were conducted. 
The first study involved examination of the response as the number of divisions 
was increased. For this purpose, was taken as 8, 16, 32, 64, 128 and 256. The 
second study involved examination of the response as the order of derivatives 
approach 1. Results of these studies are given below. 

4.1 Time Invariant FOCP 

As a first example, we consider the following Time Invariant Problem (TIP): Find 
the control u{t) which minimizes the quadratic performance index 

J{u) = l l\x\t)+u\t)]dt (19) 
2 Jo 

subjected to the system dynamics 

qD^x = -x + u, (20) 

and the initial condition 

x(0) = 1. (21) 

For this example, we have 

q(t) = r{t) = -a(t) = b{t) = xq = 1. (22) 

This example is considered here because, for a = 1, it is one of the most common 
examples of time invariant systems considered by many. The closed form solution 
for this system for a = 1 is given as (see, Agrawal (1989)) 



X 



(t) = cosh( V2t) + P sinh(y2t) (23) 



and 

u{t) = (1 + V2P) cosh(y2t) + {V2 + P) sinh(y2t) (24) 

where 

;3 = -7'''^) ^-0.9799. (25) 

y2cosh(y2) + sinh(V2) 

From Eqs. (EH) and (ESD, we get u{0) = -0.3858. 

Figures 1 and 2 show the state x{t) and the control u{t) as functions of t for 
a = 0.75 and different values of A^. 
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Figure 1: Convergence of x{t) for the TIP for a = 0.75 (A : = 8, O : = 16, 
+ : N ^32, X : N ^64:,V : N ^ 128, ★ : AT = 256) 





Or 




-0.05 












-0.1 


o 








c 
o 




o 










-0.25^ 









0.4 0.6 

Time t (sec.) 

Figure 2: Convergence of u{t) for the TIP for a = 0.75 (A : = 8, O : = 16, 
+ : AT = 32, X : AT = 64, V : A^ = 128, ★ : AT = 256) 

Figures 3 and 4 show the state x{l) and the control u{0) as a function of A^ for 
different a. From these figures, it can be seen that the solutions converge as N 
is increased, however, the convergence is slow. Further, the convergence becomes 
poor as a is decreased. Further error analysis may be necessary to identify the 
reasons for this behavior. 
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Figure 3: Convergence of a;(l) for the TIP for different a (A : a = 0.5, 
O : a = 0.75, + : a = 0.95, X:a^ 1.0) 
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Figure 4: Convergence of m(0) for the TIP for different a (A : o; 
O : q; = 0.75, + : a = 0.95, X : a ^ 1.0) 
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Figures 5 and 6 show the state x{t) and the control u{t) as functions of t for 
different values of a. These figures also show analytical results for the state x{t) 
and the control u{t) for a = 1. It can be observed that for a = 1 the numerical 
solution agrees with the analytical solution. Thus, as a approaches to 1, the 
solution for the integer order system is recovered. 
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Figure 5: State x{t) as a function of t for the TIP for different a 
( - :a = 0.5, — -.a = 0.75, ■■■ :a = 0.95, — - :a = 1) 
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Figure 6: Control u{t) as a function of t for the TIP for different a 
( - :a = 0.5, — -.a = 0.75, ■■■ :a = 0.95, -— :a = 1) 



4.2 Time Varying FOCP 

As a second example, we consider the following Time Varying Problem (TVP): 
Find the control u{t) which minimizes the quadratic performance index given in 
Eq. ( |T9l) . and which satisfies the system dynamics 

oD^x = tx + u. (26) 

The initial condition is x(0) = 1. For this example, we have 

q(t) = r{t) = hit) =xq = 1, a{t) = t. (27) 

It is one of the simplest examples of time varying systems, and for a = 1, it has 
been considered at several other places (see, Agrawal (1989), and the references 
there in). 
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Figures 7 and 8 show the state x{t) and the control u{t) as functions of t for 
different values of N. Figures 9 and 10 show the state a;(l) and the control ■u(O) 
as a function of for different a. As for the TIP, the solutions for the TVP also 
converge as N is increased, however, as before, the convergence is slow. This slow 
convergence for both examples clearly suggests that further improvement of the 
scheme is necessary. Figures 11 and 12 show the state x{t) and the control u{t) 
as functions of t for different values of a. 
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Figure 7: Convergence of x{t) for the TVP for a = 0.75 (A : = 8, 
O : N = 16, + : N = 32, X : N = 64, V : N = 128, ★ : iV = 256) 
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Figure 8: Convergence of u{t) for the TVP for a = 0.75 (A : A^ = 8, 
O : N = 16, + : N ^ 32, X : N = 6A, V : N ^ 128, ★ : A^ = 256) 
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Figure 9: Convergence of x{l) for the TVP for different a [A : a = 0.5, 
0:a = 0.75, + : a = 0.95, X : a = 1.0) 
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Figure 10: Convergence of u{0) for the TVP for different a {A : a — 0.5, 
0:a = 0.75, + : a = 0.95, X : a = 1.0) 
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Figure 11: State x{t) as a function of t for the TVP for different a 
( - :« = 0.5, — :a = 0.75, - :a = 0.95, 
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Figure 12: Control u{t) as a function of t for the TVP for different a 
( - :« = 0.5, — :a = 0.75, ••• :a = 0.95, — - :a = 1) 

This problem for o; = 1 has been solved in Agrawal (1989) using a different 
scheme.The scheme is based on the approximation with weighing coefficients and 
the lagrange multiplier technique for a class of optimal control problems (Agrawal 
(1989)). Results show that for a = 1 the numerical solutions obtained using the 
scheme developed here and in Agrawal (1989) agree well. Thus, as before, as a 
approaches to 1, the solution for the integer order system is recovered. 

It should be point out here that in integer order calculus central difference 
schemes have been used in many cases to develop numerically stable and efficient 
schemes. It is hoped that this research will initiate a similar effort in fractional 
calculus. 
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5 Conclusions 



For a general class of fractional optimal control problems a Hamiltonian was 
defined and a set of necessary conditions were derived. A direct numerical scheme 
was presented for solution of the problems. The scheme was used to solve two 
problems, time invariant and time varying. Results showed that as the number 
of divisions of the time domain was increased, the solutions converged. However, 
the convergence appears to be slow. As the value of a approaches 1, the solution 
for the integer order system is recovered. It is hoped that this research would 
initiate further research in the field, and more efficient and stable schemes would 
be found. 
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